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3 NON-LINEAR AXI SYMMETRIC POTENTIAL FLOW BOUNDARY MODEL FOR 

4 PARTIALLY CAVITATING HIGH SPEED BODIES 
5 

6 STATEMENT OF GOVERNMENT INTEREST 

7 The invention described herein may be manufactured and used 

8 by or for the Government of the United States of America for 

9 governmental purposes without the payment of any royalties 
fiO thereon or therefore . 
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J2 CROSS-REFERENCE TO RELATED PATENT APPLICATIONS 

H| 3 Not applicable. 
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3t5 BACKGROUND OF THE INVENTION 

x|'6 (1) Field of the Invention 

pk.7 The present invention relates to computer model of 

18 hydrodynamic flows and more particularly, relates to modeling 

19 partially cavitating flows over a supercavitating axisymmetric 

20 body. 

21 (2) Description of the Prior Art 

22 Modeling of boundary flows about objects subject to laminar 

23 and turbulent flows is well known in the art. High speed 

24 underwater vehicles, however, cause cavitation of the surrounding 

25 fluid. Cavitation reduces pressure in the fluid below its vapor 



1 pressure causing the fluid to vaporize, allowing the undersea 

2 vehicle to travel with lower friction when the vehicle is 

3 completely surrounded by the cavity. 

4 Partial cavitation is an unsteady phenomenon that occurs 

5 when part of the supercavitating vehicle is traveling in the 

6 cavity. Specifically, this phenomenon occurs during launch of 

7 the vehicle. A steady, partial cavitation allows development of 

8 vehicle designs which take advantage of drag reduction through 

9 cavitation. It may also be possible to take advantage of drag 
010 reduction with partial cavitation by properly directing the re- 
tf-1 entrant jet that forms in the cavity closure region. Partial 

#2 cavitation often occurs during maneuvering of the supercavitating 

yi3 vehicle. 

114 A slender body theory has been developed to solve 

JffS axisymmetric supercavitating flows. Using the slender body 

2|6 method, sources are defined along the body-cavity axis and 

T7 control points along the body-cavity surface. A nonlinear 

18 differential equation is formed by imposing dynamic boundary 

19 conditions on the cavity. A conical cavity closure is assumed in 

20 order to solve the developed nonlinear differential equation. 

21 A non-linear boundary element method for determining a 

22 cavity shape has been developed. Source and dipole strengths 

23 along the body-cavity surface are determined using kinematic 

24 boundary conditions on the wetted body surface and dynamic 

25 boundary conditions on the assumed cavity shape. The kinematic 
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1 boundary condition is then used to update the cavity shape. The 

2 process is then iterated to solve for the unknown cavity shape. 

3 Two numerical hydrodynamics models have been developed by 

4 the Naval Undersea Warfare Center for axisymmetric super 

5 cavitating high speed bodies. These models are Lhe slender body 

6 theory (SBT) model and the boundary element (BE) model. Both of 

7 these models have been proven to predict cavity shape and 

8 parameters with good accuracy. 

9 These models, however, do not account for the transition 
AO case when the vehicle is subjected to only partial cavitation. 
||1 In the SBT model, total drag is predicted by adding the 
j§2 pressure drag obtained from the model solution and the viscous 
Qp drag obtained by applying the Thwaites and Falkner-Skan 

s 14 approximations along the wetted portions of the cavitator. This 

015 method is extended to subsonic compressible flows using the 

-5 6 compressible Green's function. In the BE model, sources and 

547 dipoles are defined on the body-cavity shape and are solved using 

18- Green's formula. This yields a Fredholm integral equation of the 

19 second kind which gives the supercavitating cavity shape. 

20 Partial cavitation modeling has been done by Uhlman, J.S. 

21 (1987), The Surface Singularity Method Applied to Partially 

22 Cavitating Hydrofoils, Journal of Ship Research, Vol. 31, No. 2, 

23 pp. 107-24; Uhlman, J.S. (1989), The Surface Singularity or 

24 Boundary Integral Method Applied to Supercavitating Hydrofoils, 

25 Journal of Ship Research, Vol. 33, No. 1, pp. 16-20; Kinnas, 
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1 S.A., and Fine, N.E. (1990), Non-Linear Analysis of the Flow 

2 Around Partially and Super-Cavitating Hydrofoils by a Potential 

3 Based Panel Method, Proceedings of the IABEM-90 Symposium, 

4 International Association for Boundary Element Methods, Rome, 

5 Italy, and Kinnas, S.A., and Fine, N.E. (1993), A Numerical 

6 Nonlinear Analysis of the Flow Around Two- and Three-Dimensional 

7 Partially Cavitating Hydrofoils, Journal of Fluid Mechanics, Vol. 

8 254. However, these methods are explicitly adapted for 

9 hydrofoils, and the theories presented therein are not readily 
QfO adapted to supercavitating vehicles. 

3J1 

•152 SUMMARY OF THE INVENTION 

Q3 One object of the present invention is a method for modeling 

^14 partial cavitation. 

E|5 Another object is that such method model partial cavitation 

j|6 about an axisymmetric vehicle. 

H>7 Accordingly, the present invention provides a method for 

18 calculating cavity shape for partial cavities about an 

19 axisymmetric body having a cavitator located at the foremost end. 

20 The method includes receiving system parameter data including 

21 geometric data describing the axisymmetric body, a cavity length, 

22 and a convergence tolerance. Boundary element panels are 

23 distributed along the body-cavity surface and matrices are 

24 initialized for each boundary element panel using the unit 

25 dipole, unit source functions and known boundary values. 
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1 Disturbance potential matrices are formulated for each boundary 

2 element panel using disturbance potentials, normal derivatives of 

3 disturbance potentials, and no net flux boundary conditions. The 

4 initialized matrices and the formulated matrices are solved for 

5 each boundary panel to obtain unknown disturbance potentials 

6 along the wetted body-cavity surfaces, and normal derivatives of 

7 disturbance potentials along the cavity surface. The cavity 

8 position is then updated by moving each panel to satisfy the 

9 kinematic boundary condition, no flux across the cavity. The ^ 
^tfo method then tests for convergence against a tolerance, and steps 
y|l are iterated until convergence is achieved. The method then 

Jl2 provides parameters of interest and the location of the cavity as 

output. Another aspect of this invention allows the calculation 

%4 of cavity shape and cavity length for an input cavitation number. 

35 This is accomplished by an outer loop adjusting cavity length 

F |6 until the model converges to the input cavitation number. 
^17 

18 BRIEF DESCRIPTION OF THE DRAWINGS 

19 These and other features and advantages of the present 

20 invention will be better understood in view of the following 

21 description of the invention taken together with the drawings 

22 wherein: 

23 FIG. 1 is a diagram of a partially cavitating axisymmetric 

24 body related to the method of the current invention; and 
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1 FIG. 2 is a flow chart of the method of the current 

2 invention; and 

3 FIG. 3 is a flow chart of another method of the current 

4 invention, 

5 

6 DESCRIPTION OF THE PREFERRED EMBODIMENT 

7 FIG. 1 shows a diagram of the physical problem of partial 

8 cavitation. FIG. 1 shows a radial cross section of an 

9 axisymmetric body 10. Axis r represents the radius from the axis 
r3j0 of body 10. Axis x represents the length along the body 10 

fjjil measured from a cavitator disk 12. Although a cavitator disk is 

4E2 shown, the model can calculate cavities for cavitator cones as 

Q3 well as cavitator disks. Flow, U w , is in the direction of arrow 

!14 14. A cavity 16 is shown extending from the edge of the 

2"5 cavitator along the length of body 10. The length of the cavity, 

.fl6 l cf is shown by dimension arrows. Likewise, the length of the 

17 body, l hf is also shown by dimension arrows. 

18 Body 10 extends beyond a cavity closure 18. Cavity 16 is 

19 closed to the body 10 with a modified Riabouchinsky cavity 

20 termination wall. Cavity closure 18 can be positioned in either 

21 body conical section 22 or body cylindrical section 24. The 

22 plane of cavity closure 18 is referenced in the following 

23 disclosure as an endplate. 

24 Body 10 has a flat front area 20 followed by a conical 

25 section 22 and a cylindrical section 24. The diameter of flat 
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1 front area 20 should be less than or equal to the diameter of the 

2 cavitator disk 12 base. 

3 The mathematical formulations in of this algorithm are based 

4 on using the cavitator diameter to remove dimensionality for all 

5 lengths and using the free stream velocity, U w , to remove 

6 dimensionality for all velocities. Alternate formulations using 

7 standard units can also be developed. 

8 The flow field is governed by Laplace's equation, 

9 V 2 <D = 0 (1) 

Qlo where $ is the total potential which is the sum of free 

3jl stream potential, (j^ , and disturbance potential, <J>, giving: 

U2 3> = (L+4) (2) 

3|3 The free stream potential is the product of the velocity and 

{34 the distance, x. Because the equation has been non- 

r$5 dimensionalized, the velocity is 1, and the free stream 

T16 potential, (j)^ , is x. The disturbance potential, <j>, also obeys 

17 Laplace's equation, giving: 

18 V 2 ^ = 0 (3) 

19 The disturbance potential satisfies Green's third identity, 

20 yielding a Fredholm integral equation of the second kind along 

21 the cavitator, cavity, endplate and body. Thus, at any point, x, 

22 on the body-cavity surface, the disturbance potential can be 
2 3 computed from: 
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27T<t>(x) = 



(|)(x) — G(x; x') (|)(x)G(x; x') 

8n dn 



dS (4) 



where x ! are the points where the sources and dipoles are 
distributed under the boundary element model; 
S is the body-cavity surface; and 
G(x,x T ) is the Green function. 
The Green function is further identified as: 

G(x,x f ) = —^— (5) 
\x — x I 

The dynamic condition on the cavity boundary is derived from 
Bernoulli's equation. Along the cavity surface, this can be 
written as: 

p x +j P Ui=p c +±pU 2 s (6) 

where jt^is the free stream ambient pressure; 

p is the free field fluid density; 

p c is the pressure inside the cavity; and 

U s is the flow velocity at the cavity surface. 

The flow velocity at the cavity surface can be obtained from 

equation (6) giving: 



U s = VT+a ( 7 ) 

where a is the cavitation number which is defined as: 

<T = ^^ (8) 

The kinetic boundary condition is that no flow crosses the body- 



1 cavity boundary, 

on 

3 where n x is the axisymmetric body free-stream velocity power. 

4 The no net flux condition, 

5 H°MdS = 0 (10) 

» dn 

6 is also required to make the problem a determinate system.. 

7 Total drag is calculated by adding the drag coefficients. 
_8 The pressure drag coefficient, C p , at x is calculated as 

5f9 follows: 

|0 C p = l-U(x) 2 (11) 

r§l The pressure contribution to the drag coefficient may then be 

s 12 computed as : 

S3 C dp =-$C p n x dS (12) 

ri4 The viscous contribution to the drag coefficient along the wetted 

15 portions of the conical and cylindrical body areas is calculated 

16 using the International Towing Tank Conference equation given by 

17 Newman, Marine Hydrodynamics , MIT Press, Cambridge, Mass. 1980, 

18 for the friction coefficient, C f , at x is as follows: 

0.075 
f (log 10 (7?(x)-2)) 

2 0 where R(x) is the local Reynolds number. 

21 The total viscous drag coefficient, C dv , is: 
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1 C dv =-Hc f s x dS (14) 

2 The base drag coefficient, C db , which is the component of 

3 pressure drag associated with the base of the body is: 

5 where b base is the body radius at the base. 

6 The total drag coefficient is then given by 

7 C,=C++C„+C db . (16) 

^8 The panels are distributed along the cavitator, cavity,, endplate, 

2^9 and cylindrical body section aft of the cavity, according to the 

r^O partial floor method, known in the art. The partial floor method 

3.1 optimizes the number of panels in accordance with requirements 

pi2 for getting good convergence. Non-uniform panel spacing is used 

33 in many locations, in order to reduce the number of panels 

f§4 without reducing the accuracy of the solution. 

15 During iteration, the end plate height is determined by 

16 integrating the cavity surface back from its detachment point on 

17 the cavitator, and the number and distribution of panels along 

18 the endplate changes according to the changes in the endplate 

19 height. Smaller panels are required at highly non-linear flow 

20 locations, such as the region near the cavitator. Panel 

21 distribution in the wetted body area after cavity closure 18 

22 changes to keep the aspect ratio of the neighboring panels 
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1 between 0.5 and 2.0, in order to ensure good accuracy of the 

2 results. 

3 In following the method of the current invention, first an 

4 initial cavity is defined. An arbitrary initial cavity can be 

5 chosen as a cone extending from the cavitator edge to an assumed 

6 endplate height of 0.2 or 0.3 is sufficient for most cases. In 

7 this discussion, the endplate height is measured as the radial 

8 offset from the body surface to the last point of the cavity. By 

9 applying equation (4) on all panels along the cavity body 

t§0 surface, S, a system of equations is obtained. This system is 

^{1 solved for the disturbance potentials, ((>, along the wetted 

[52 portions of the boundary and on the Riabouchinsky endplate; the 

33 normal derivative of the disturbance potential along the cavity 

fl;4 boundary; and the cavitation number. 

§5 The kinetic boundary condition given in equation (9) is 

E§6- applied along cavitator, endplate, and aft body to update the 

jiass 

'17 cavity shape. In order to update the cavity, the program 

18 calculates how much each panel has to be rotated to satisfy the 

19 no flow condition. The program starts with the first panel at 

20 the cavitator and shifts the aft most point of the panel in the 

21 radial direction which satisfies the calculated rotation. The 

22 panel is rotated with the aft most point. The foremost point of 

23 the next panel is then shifted to the same radius as the previous 

24 aft most point. This process is continued until the panel 

25 adjacent to the endplate is updated. The endplate height is 
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1 adjusted to the aft most point of the aft cavity panel. The 

2 iteration continues until the kinetic boundary condition 

3 converges to within a tolerance, giving the cavity shape. 

4 From the converged disturbance potential along S, the 

5 disturbance velocity components can be calculated: 

6 u x =^ and u r =^. (17) 

8x dr 

7 Referring now to FIG. 2, there is shown a flowchart of the 

8 current invention. In the input step 30, geometric and other 

pi 9 system parameter data including the estimated cavitation number, 

||0 the estimated cavity length and the convergence criteria is read, 

jfl The routine then distributes boundary element panels along the 

£|2 cavitator, cavity, endplate, body extension in the conical 

s 13 section, body extension in the horizontal section, and the aft 

B|4 body. The panels are distributed in order to reduce the number 

S|5 of panels and get an accurate result. In the initialize step 32, 

6 the algorithm calculates the unit dipole and unit source 

17 functions and initializes matrices for the influence functions 

18 with known boundary values wherever applicable. The formulate 

19 equations step 34 formulates matrices for each panel using the 

20 disturbance potential equation (4) and no net flux condition 

21 given in equation (9) . The solve equations step 36 solves the 

22 matrices created in the formulate equations step 34 in order to 

23 obtain the unknown disturbance potential along wetted body 

24 sources, normal distributions of disturbance potentials along 

25 cavity surfaces, and the cavitation number. The compute forces 
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1 step 38 computes velocity components such as those in equation 

2 (17) and drag coefficients: including pressure drag, equation 

3 (12); viscous drag, equation (14); and base drag, equation (15) 

4 from the solved equations. In the update cavity step 40, the 

5 cavity is updated from the computed forces using the kinetic 

6 boundary condition of equation (9). Convergence on cavity shape 

7 is checked in the converges decision step 42. If the cavity is 

8 not converged, the initialize step 32 is executed to calculate 

9 influence functions for the updated cavity and next iteration 
Q|0 thus begins. Once the cavity has converged, the compute 

fill parameters step 4 4 computes various output parameters of the 

M2 converged solution which include pressure drag, viscous drag, 

33 base drag, total drag, cavitation number, cavity length, maximum 

-14 cavity radius, length of cavity to maximum radius location. The 

E§5 output results step 46 then provides the location of the cavity 

116 written as coordinates and the cavity's disturbance potential, 

Ht7 disturbance potential gradient, and pressure coefficient. 

18 The basic algorithm enumerated above provides cavity shape 

19 and cavitation number based on an input cavity length. In order 

20 to obtain cavity shape and cavity length for an input cavitation 

21 number, the embodiment of FIG. 3 adds an additional series of 

22 iterations. The user inputs a cavitation number and an assumed 

23 cavity length. This embodiment follows the previous embodiment 

24 in converging on a new cavitation number, a, for the assumed 

25 cavity length. In step 48, if the new cavitation number is 
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1 within a tolerance of the given cavitation number, parameters are 

2 computed, step 44, and the results are provided, step 46* 

3 Otherwise the embodiment proceeds to step 50 wherein the 

4 algorithm determines the relationship between the new cavitation 

5 number, a, and the given cavitation number. In step 52, cavity 

6 length is increased by a predetermined amount if the calculated 

7 cavitation number is lower than the initial cavitation number, 

8 and in step 54 the cavity length is decreased by a predetermined 

9 amount if the calculated cavitation number is greater than the 
30 initial cavitation number. The routine loops back to the 

Offl initialize step 32 and recalculates the cavitation number for the 

4S2 new cavity length. Operation continues until the calculated 

Cl3 cavitation number falls within a tolerance of the initial 

=14 cavitation number, the cavity length has converged, as tested in 

J5 step 48. 

r j6 Using this invention, partial cavitation for high-speed 

"HE 7 underwater bodies can be analyz'ed. As disclosed, the invention 

18 can analyze axisymmetrical bodies using two cavitator shapes, a 

19 disk and a cone; however, the invention can easily be modified to 

20 analyze other axisymmetric cavitator shapes. As disclosed the 

21 inventive method can converge on cavity length or cavitation 

22 number. Total drag is calculated by adding the pressure drag, 

23 viscous drag and base drag. The invention can also be utilized 

24 for studying the effects of body aft radius, body cone angle and 

25 body cone angle starting at the cavity closure if the closure is 
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1 on conical section 22. This method provides new information 

2 concerning the physics of cavitation which can be used in the 

3 . design of cavitating vehicles. 

4 In light of the above, it is therefore understood that 

5 within the scope of the appended claims, the invention may be 

6 practiced otherwise than as specifically described. 
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What is claimed is: 



1. A method for calculating parameters for an axisymmetric 
partially cavitating body having a cavitator located at the 
foremost end, said method comprising the steps of: 

receiving system parameter data including geometric data 
describing the axisymmetric body, a convergence 
tolerance, and an initial cavity shape including an 
endplate height, an endplate location, and a cavity 
length; 

initially distributing boundary element panels along the 
initial cavity shape, endplate and the axisymmetric 
body aft of the endplate; 

initializing matrices for each boundary element panel using 
the disturbance potentials at the boundary element 
panels and known boundary values; 

formulating disturbance potential matrices for each boundary 
element panel utilizing disturbance potential equations 
and no net flux boundary conditions; 

solving initialized matrices and formulated disturbance 

potential matrices for each boundary panel to obtain 
unknown panel sources, unknown dipoles and unknown 
cavitation numbers ; 
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computing forces and velocities at each panel from the panel 
sources, dipoles and cavitation numbers to obtain 
velocity components, pressure drag, viscous drag, and 
base drag; 

updating the cavity by moving each panel in accordance with 
the calculated forces and velocities and the boundary 
conditions; 

testing for convergence by comparing the movement of each 
panel against the convergence tolerance; 

iterating said steps of initializing matrices, formulating 

matrices, solving formulated matrices, computing forces 
and velocities, updating the cavity and testing for 
convergence when said test for convergence indicates 
that movement of at: least one panel exceeds the 
convergence tolerance; 

computing parameters of interest when said test for 

convergence indicates that movement of all panels is 
within the convergence tolerance; and 

outputting the location of the cavity and the computed 
parameters . 
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